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Abstract 



A discrete time control algorithm using the damped least squares is 
introduced for acceleration and energy exchange controls in nonlinear vi- 
brating systems. It is shown that the damping constant of least squares 
and sampling time step of the controller must be inversely related to insure 
that vanishing the time step has little effect on the results. The algorithm 
is illustrated on two linearly coupled DufHng oscillators near the 1:1 in- 
ternal resonance. In particular, it is shown that varying the dissipation 
ratio of one of the two oscillators can significantly suppress the nonlinear 
beat phenomenon. 

1 Introduction 

The damped least squares is a simple but effective analytical manipulation that 
helps to avoid singularity in practical minimization and control algorithms. It 
is also knovifn as Levenberg-Marquardt method [llj. In order to illustrate the 
idea in simple terms, let us consider the minimization problem 



where E £ i?" is a given vector, the notation ||...|| indicates the Euclidean 
norm in i?", A is typically a Jacobian matrix of n rows and m columns, and 
6u e i?™ is an unknown minimization vector. Although a formal solution of this 
problem is given by Su — {A'^A)^^A'^E, the matrix product A^A may appear 
to be singular so that no unique solution is possible. This fact usually points 
to multiple possibilities of achieving the same result unless specific conditions 
are imposed on the vector Su. The idea of damped least squares is to avoid 
such conditioning by adding one more quadratic form to the left hand side of 
expression (IT]) as follows 



\\E — yl(57i||^ — > min 



(1) 



\\E - ASu\\^ + X\\Suf min 



(2) 



where A is a positive scalar number, which is often called damping constant] 
note that the term 'damping' has no relation to the physical damping or energy 
dissipation effects in vibrating systems usually characterized by damping ratios. 

Now the inverse matrix includes the damping constant A which can provide 
the uniqueness of solution given by 



where / is n x n identity matrix. 

Different arguments are discussed in the literature regarding the use of 
damped least squares and best choice for the damping parameter A [T], 0, 
[3], [4], [6], [7], 0, [lO], [H], [16], [n], [23], [24]. In particular, it was noticed 
that the parameter A may affect convergence properties of the corresponding 
algorithms. The parameter A can be used also for other reason such as shifting 
the solution 5u into desired area in i?™. In this case, the meaning of A is rather 
close to that of Lagrangian multiplier imposing constraints on control inputs. 

In case of dynamical systems, when all the quantities in ^ may depend on 
time, a continuous time analogue of ([2|) can be written in the integral form 



where the interval of integration is manipulated as needed, for instance, T can 
be equal to sampling time of the controller [12j . 

However, in the present work, a discrete time algorithm based on the damped 
least squares solution ([s]) , which is used locally at every sample time i„ , is intro- 
duced. Such algorithm appears to be essentially discrete namely using different 
time step h may lead to different results. Nevertheless, if the parameters A and 
h are coupled by some condition then the control input and system response 
show no significant dependence on the time step. 

A motivation for the present work is as follows. In order to comply with 
the standard tool of dynamical systems dealing with differential equations, the 
methods of control are often formulated in continuous time by silently assuming 
that a discrete time analogous is easy to obtain one way or another whenever it is 
needed for practical reasons. For instance, data acquisition cards and on-board 
computers of ground vehicles usually acquire and process data once per 0.01 sec. 
Typically, based on the information, which is known about the system dynamic 
states and control inputs by the time instance the computer must calculate 
control adjustments for the next active time instance, i„+i. The corresponding 
computational time should not therefore exceed t„+i — tn — 0.01 sec. Generally 
speaking, it is possible to memorize snapshots of the dynamic states and control 
inputs at some of the previous times {..., t„_2, However, increasing the 

volume of input data may complicate the code and, as a result, slow down the 
calculation process. Therefore, let us assume that updates for the control inputs 
are obtained by processing the system states, controls, and target states given 
only at the current time instance, tn- The corresponding algorithm can be built 
on the system model described by its differential equations of motion and some 



5u = [A 



^A + \I)-^A^E 



(3) 
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rule for minimizing the deviation (error) of the current dynamic states from the 
target. Recall that, in the present work, such a rule will be defined according 
to the damped least squares Illustrating physical example of two linearly 
coupled Duffing oscillators is considered. It is shown that the corresponding 
algorithm, which is naturally designed and effectively working in discrete time, 
may face a problem of transition to the continuous time limit. 



2 Problem formulation 

Consider the dynamical system 

x^f{x,x,t,u) (5) 

where x = x{t) G i?" is the system position (configuration) vector, the overdot 
indicates derivative with respect to time t, the right-hand side / S i?" represents 
a vector-function that may be interpreted as a force per unit mass of the system, 
and u = u{t) e i?™ is a control vector, whose dimension may differ from that 
of the positional vector x so that generally n ^ m. 

In common words, the purpose of control u{t) is to keep the acceleration 
x{t) of system ([s]) as close as possible to the target x*{t). The term 'close' 
will be interpreted below through a specifically designed target function of the 
following error vector 

E[t) = i*{t)~-i(t) (6) 

As discussed in Introduction, for practical implementations, the problem 
must be formulated in terms of the discrete time {tk} as follows. Let Xk = 
x{tk), Xk = i(ifc), and = u(tk) are observed at some time instance t^. The 
corresponding target acceleration, = x*(tk), is assumed to be known. Then, 
taking into account ^ and Q, gives the following error at the same time 
instance 

Ek ^ xl - f{xk,Xk,tk,Uk) (7) 
Now the purpose of control is to minimize the following target function 

Pk = \ElWkEk (8) 

= ^[aifc - fixk.Xk.tk.Ukyi^Wklxl - f{xk,Xk,tk,Uk)] 

where Wk is n x n diagonal weight matrix whose elements are positive or at 
least non-negative functions of the system states, Wk — W{xk,Xk,tk)- 

Note that all the quantities in expression ([s]) represent a snapshot of the 
system ai t = tk while including no data from the previous time step tk^i- 
Since the control vector Uk cannot be already changed at time tk then quantity 
Pk is out of control at time tk- In other words expression ([s]) summarizes all 
what is observed now, at the time instance tk- The question is how to adjust 
the control vector u for the next step tk+i based on the information included in 
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([8| while the system state at t — tk+i is yet unknown, and no information from 
the previous times {..., i„_2, tn-i} is available. 

Let us represent such an update for the control vector in the form 

Uk+i = Uk + Suk (9) 

were Suk is an unknown adjustment of the control input. 
Replacing Uk in ([s]) by (|9]) and taking into account that 

f{xk,ik,tk,Uk+i) = f{xk,ik,tk,Uk)+Ak6uk+0{\\Suk\\^) (10) 
Ak = df{xk,ik,tk,Uk)/duk 



gives 



Pk ^\{Ek-' Ak5ukfWk{Ek - Ak6uk) 



(11) 



where Ak is the Jacobian matrix of n rows and m columns. 

Although the replacement Uk by Uk+i in ( 10 1 may look artificial, this is how 
the update rule for the control vector u is actually defined here. Namely, if 
Uk did not provide a minimum for Pk{xl.,Xk,Xk,tk,Uk), then let us minimize 
Pk{x1, Xk,Xk,tk,Uk + Suk) with respect to 5uk and then apply the adjusted 
vector ([9]) at least the next next time, tn+i- Assuming that the variation Suk is 



small, in other words, Uk is still close enough to the minimum, expansion ( 10 1 



is applied. Now the problem is formulated as a minimization of the quadratic 
form (11) with respect to the adjustment Suk- However, what often happens 



practically is that function (111 has no unique minimum so that equation 



dPk 
d6uk 



= 



(12) 



has no unique solution. In addition, even if the unique solution does exist, it 
may not satisfy some conditions imposed on the control input due to the physical 
specifics of actuators. As a result, some constraint conditions may appear to 
be necessary to impose on the variation of control adjustment, 6uk- However, 
the presence of constraints would drastically complicate the problem. Instead, 



the target function (11) can be modified in order to move solution Suk into the 



allowed domain. For that reason, let us generalize function (11) as 



-(Ek - AkSukfWkiEk - AkSuk) 
+ l{Bk + CkSukfAkiBk + CkSuk) 



(13) 



where A^ = A{xk, Xk,tk) is a diagonal regularization matrix, Bk — B{xk, Xk,tk) 
is a vector-function of n elements, and Ck ~ C{xk,Xk, tk) is a matrix of n rows 
and m columns. 

Note that the structure of new function (13) is a generalization of ([2]). Sub- 
stituting (13) in (12), gives a linear set of equations in the matrix form whose 
solution 5uk brings relationship ^ to the form 
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Uk+1 =uk + {AlWkAk + ClKkCk)-\AlWkEk - C^AkEk) (14) 



The entire discrete time system is obtained by adding a discrete version 
of the dynamical system ^ to ( 14 ) . Assuming that the time step is fixed, 
tk+i — tk — h, a simple discrete version can be obtained by means of Euler 
explicit scheme as follows 



Vk+i = Vk + hf{xk,Vk,tk,Uk) 



(15) 



Finally, equations ( 14 1 and ( 15 ) represent a discrete time dynamical system, 
whose motion should follow the target acceleration xl = x*{tk). 



It will be shown in the next section that the structure of equation ( 14 ) does 



not allow for the transition to continuous limit of the entire dynamic system ( 14 1 



through (15), unless some specific assumption are imposed on the parameters 
in order to guarantee that 5uk — 0{h) as — > 0. 



3 The illustrating example 

The algorithm, which is designed in the previous section, is applied now to a 
two-degrees-of-freedom nonlinear vibrating system for an active control of the 
energy exchange (nonlinear beats) between the two oscillators. The problem of 
passive control of energy flows in vibrating systems is of great interest [22] , and 
it is actively discussed from the standpoint of nonlinear beat phenomena |14j . 
The beating phenomenon takes place when frequencies of the corresponding 
linear oscillators are either equal or at least close enough to each other. 

For illustrating purposes, let us consider two unit-mass Duffing oscillators of 
the same linear stiffness K coupled by the linear spring of stiffness 7. The system 
position is described by the vector- function of coordinates, x[t) = {xi{t), X2{t))'^ . 
Introducing the parameters 17 = (7 -|- Ky/"^ and £ = 7/(7 + K), brings the 
differential equations of motion to the form 

Xi = Vl 
X2 = V2 

ii = -2Cnvi-n'^xi+e{n'^X2~axl) = fi{xi,X2,vi) (16) 

V2 = —2uil,V2 — ^'^X2+s{V,'^Xi~aX2)^f2{xi,X2,V2,U) 

where a is a positive parameter, ( and u are damping ratiot[^of the first and the 
second oscillators, respectively; the damping ratio u, which is explicitly shown 
as an argument of the function f2ixi,X2,V2,u), will be considered as a control 
input. 

^As mentioned in Introduction, the damping (dissipation) ratio should not be confused 
with the damping coefficient A. 
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The problem now is to find such variable damping ratio u — u{t) under 
which the second oscillator accelerates as close as possible to the given (target) 
acceleration, X2{t). 

Following the discussion of the previous section, let us consider the prob- 
lem in the discrete time {tk}- In order to avoid confusion, the iterator k 
will be separated from the vector component indexes by coma, for instance, 
Xk = (xi^k,X2,k)'^ ■ Since only the second mass acceleration is of interest and 
the system under consideration includes only one control input w, then, assum- 
ing the weights to be constant, gives 



Wk = 




1 



Ak = 



d 
duk 



fl.k 

f2,k 



where fi,k = fi{xi^k,X2,k,vi,k) and /2,/c = f2ixi^k,X2,k,V2,k,Uk), and other 
matrix terms become scalar quantities, say, A^; — A, Bk — b, and Ck = 1- The 
unities in Wk and Ck can always be achieved by re-scaling the target function 
and parameters A and b. Note that re-scaling the target function by a constant 



factor has no effect on the solution of equation ( 12 ) 



As a result, the target function (13) takes the form 



Pk = 



• * r 9f2.k r 

^9 h — 19 k 7^ — —OUk 

2.k J2.k k 



Suk? 



(17) 



In this case, equation ( 12 ) represents a single linear equation with respect 
to the scalar control adjustment, Suk- Substituting the corresponding solution 
in (14 1 and taking into account (|15|), gives the discrete time dynamical system 



Mfc+l = Uk 



{.f2.k-xlk)idf2.k/duk) 
{dh.k/d^k? + A 



\b 



(18) 



and 



Xl.k+l 
X2,k+l 
Vl.k+1 
V2,k+1 



xi.k + hvi^k 

X2.k + hV2,k 

vi.k + hfi^k 

V2.k + hj2k 



(19) 



Let us assume now that the target acceleration X2 is zero, in other words, 
the purpose of control is to minimize acceleration of the second oscillator at any 
sample time tk as much as possible. Let us set still arbitrary parameter b also 



to zero. Then the target function ( 17 1 and dynamical system ( 18 ) and ( 19 1 take 
the form 



Pk = 



1 



f2ixi,k,X2,k,V2,k,Uk) + 



df2ixi^k,X2,k,V2,k,Uk] 

duk 



Suk 



+ ^i6ukr (20) 
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Uk+1 = Uk 



2VLV2.. 



Xl,k+1 
X2,k+1 
Vl,k+1 
V2,k+1 



xi^k + hvi^k 

X2,k + hv2^k 

Vl.k + hfi{xi^k,X2,k,Vi^k) 
V2,k + hf2{xi^k,X2,k,V2,k,Uk) 



(21) 



where the functions /i and /2 are defined in ( 16 1 



As fohows from the first equation in (21), transition to the continuous time 



hmit for the entire system (|2l| would be possible under the condition that 

= 0{h), 



2Jlt>2 t 



as 







(22) 



Condition ([22| can be satisfied by assuming that fl — 0{h). Such an as- 
sumption, however, makes little if any physical sense. As an alternative choice, 
the condition A = 0{h^^) can be imposed by setting, for instance, 



Xh — Ao 



(23) 



where Aq remains finite as — ?> 0. 



However, condition (231 essentially shifts the weight on control to the second 



term of the target function (17) so that the function asymptotically takes the 
form 



Pk 



Ao 
2h 



{5ukf 







(24) 



Such a target function leads to the solution 5uk — 0, which effectively elim- 
inates the control equation. In other words, the iterative algorithm seems to 
be essentially discrete. As a result, the control input Uk, generated by the first 



equation in (21), depends upon sampling time interval h. Let us illustrate this 



observation by implementing the iterations (21 ) under the fixed set of parame- 
ters, £ = 0.1, ri = 1.0, a — 1.5, C — 0.025, and initial conditions, uq — 0.025, 
xifi = 1.0, X2fi = 0.1, vifl = V2.0 = 0. The values to vary are two different sam- 
pling time intervals, h = 0.01 and h = 0.001, and three different values of the 
damping constant, A = 0.1, A = 1.0, and A = 10.0. For comparison reason. Fig. 
1 shows time histories of the system coordinates under the fixed control vari- 
able u — This (no control) case corresponds to free vibrations of the model 



( 16 ) whose dynamics represent a typical beat- wise decaying energy exchange 
between the two oscillators. As mentioned at the beginning of this section, the 
beats are due to the 1:1 resonance in the generating system {e ^ 0, u ~ ( — 0); 
more details on non-linear features of this phenomenon, the related analytical 
tools, and literature overview can be found in [2D] and |14| . In particular, the 



standard averaging method was applied to the no damping case of system ( 16 ) 
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Now the problem is to suppress the beat phenomenon by preventing the 
energy flow from the first osciUator into the second osciUator. As follows from 
Figs. 2 through 5, such a goal can be achieved by varying the damping ratio 
of the second oscillator, {uk}, during the vibration process according to the 
algorithrrj^ (21). First, the diagrams in Figs. 2 and 3 confirm that the sampling 



time interval h represents an essential parameter of the entire control loop. In 
particular, decreasing the sampling interval from h = 0.01 to h = 0.001 effec- 
tively increases the strength of the control; compare fragments (b) in Figs. 2 
and 3. However, if such decrease of the sampling time is accompanied by the 



increase of A according to condition (23), then the strength of control remains 



practically unchanged; compare now fragments (b) in Figs. 2 and 4. As follows 
from fragments (a) in Figs. 2 and 4, the above modification of both parameters, 
h and A, also brings some difference in the system response during the interval 
80 <t < 150, but this is rather due to numerical effect of the time step. Finally, 
analyzing the diagrams in Figs. 3 and 5, shows that reducing the parameter A 
as many as ten times under the fixed time step h leads to a significant increase 
of the control input {uk} with a minor effect on the system response though. 
Therefore the parameter A can be used for the purpose of satisfying some con- 
straint conditions on the control inputs {uk} in case such conditions are due 
to physical limits of the corresponding actuators. In addition, let us show that 



parameter A may affect the convergence of algorithm (21 ) based on the following 
convergence criterion |18j : 

For a fixed point to be a point of attraction of the algorithm Zk+i = G{zk) 
a sufficient condition is that the Jacobian matrix of G at the point has all 
its eigenvalues numerically less than 1, and a necessary condition is that they 
are numerically at most 1. The geometric rate of convergence is the numerically 
largest eigenvalue of this Jacobian. 



Applying this criterion to the algorithm (21) at zero point, gives that one 
of the eigenvalues is always zero, qo = 0, whereas another four eigenvalues, qi 
(i = 1, ...,4) are proportional to the time step, qi = hpi, where the coefficients 
Pi are given by the roots of algebraic equation 

/ -I- 2Cnp^ + 2r2V + 2Cfl^p + (1 - e2)si4 ^ (25) 



As follows from ( 25 ) , the damping coefficient A has no influence on the con- 
vergence condition near the equilibrium point, and the convergence can always 
be achieved under a small enough time step h. Nevertheless, the damping coef- 
ficient may appear to affect the convergence away from the equilibrium point. 
In this case, analytical estimates for eigen values of the Jacobian become tech- 
nically complicated unless £ = 0, when four of the five eigenvalues vanish as 
/i — >■ 0, except one eigenvalue, which is estimated by 

<i = -(^ + :A:2) (26) 



^Note that, although the algorithm is designed to suppress accelerations of the second 
oscillator, acceleration and energy levels of vibrating systems are related. 
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This root gives g — > go = as U2 0. However, when V2 0, equation 
(26) gives the estimate 0<(7<lasoo>A>0. Therefore, only the necessary 
convergence condition is satisfied for A = 0. 



4 Conclusions 

In this work, a discrete time control algorithm for nonlinear vibrating systems 
using the damped least squares is introduced. It is shown that the corresponding 
damping constant A and sampling time step h must be coupled by the condition 
Xh = constant in order to preserve the result of calculation when varying the 
time step. In particular, the above condition prohibits a direct transition to 
the continues time limit. This conclusion and other specifics of the algorithm 
are illustrated on the nonlinear two-degrees-of-freedom vibrating system in the 
neighborhood of 1:1 resonance. It is shown that the dissipation ratio of one of 
the two oscillators can be controlled in such way that prevents the energy ex- 
change (beats) between the oscillators. From practical standpoint, controlling 
the dissipation ratio can be implemented by using devices based on the physical 
properties of magnetorheological fluids (MRF) [8], [19]. In particular, differ- 
ent MRF dampers are suggested to use for semi-active ride controls of ground 
vehicles and seismic response reduction. 
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Figure 2: Beat suppression under the time increment h = 0.01 and weight 
parameter A = 1.0: (a) the system response, (b) control input - the damping 
ratio of second oscillator. 
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Figure 3: Beat suppression under the reduced time increment h = 0.001 and 
the same weight parameter A = 1.0: (a) the system response, (b) control input 
- the damping ratio of second oscillator. 
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Figure 4: Beat suppression under the reduced time increment h = 0.001 but 
increased weight parameter A = 10.0: (a) the system response, (b) control input 
- the damping ratio of second oscillator. 
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Figure 5: Beat suppression under the reduced time increment h — 0.001 and 
vanishing weight parameter A = 0.1: (a) the system response, (b) control input 
- the damping ratio of second osciUator. 
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